Electronic structure of noble metal impurities in semiconductors: Cu in GaP 
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A numerical method for calculation of the electronic structure of transition metal impurities in 
semiconductors based on the Green function technique is developed. The electronic structure of 3d 
impurity is calculated within the LDA+U version of density functional method, whereas the host 
electron Green function is calculated by using the linearized augmented plane wave expansion. The 
method is applied to the Cu impurity in GaP. The results of calculations are compared with those 
obtained within the supercell LDA procedure. It is shown that in the Green function approach Cu 
impurity has an unfilled 3d shell. This result paves a way to explanation of the magnetic order in 
dilute Gai-ajCu^P alloys. 

PACS numbers: 71.55.Eq,71.15.Ap,75.50.Pp 
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I. INTRODUCTION 

The experimental and theoretical studies of dilute 
magnetic semiconductors (see, e.g. the recent review^) 
have revived the interest to details of reconstruction of 
the electronic structure of host materials induced by tran- 
sition metal ions and concomitant defects. This inter- 
est stems from the fact that the simple Vonsovskii-Zener 
model of s-d exchange is apparently not sufficient for an 
exhaustive explanation of the behavior of the most pop- 
ular system (Ga,Mn)As, 2 not to mention the wide-gap 
materials like (Ga,Mn)N, (Zn,Co)0, (Ti,Co)0 2 Not 
only the localized spin of magnetic ions but also the ac- 
ceptor or donor-like states in the energy gap related to 
these ions are involved in the indirect exchange between 
the magnetic ions responsible for the long-range magnetic 
order. The nature of these states is the matter of a vivid 
discussion in the current literature. 

In particular, an isolated Mn impurity in GaAs cre- 
ates a 0.11 eV acceptor level relative to the top of va- 
lence band. Besides, the electrons in the half-filled 3d 
shell form resonance levels in the middle of this band 
because of an anomalous stability of the half-filled 3d 5 
shell&i Since the substitution impurity Mn 2+ (3d 5 ) is 
negatively charged relative to the host semiconductor, 
localizing a hole makes this defect neutral, and the bind- 
ing energy of this hole is provided by the combined action 
of the Coulomb potential, central cell substitution poten- 
tial, hybridization and, maybe, s-d exchange*^ At a high 
enough Mn concentration, these acceptor levels form an 
impurity band and eventually merge with the hole states 
near the top of the valence band (see Ref. [l(] for a de- 
tailed discussion of the current experimental situation). 

According to the available calculations of the electronic 
spectra of an isolated Cu in GaPfii the copper impurity 
should have a similar electronic structure. Due to the 
special stability of the filled 3d 10 shell all the 3d levels 
of the Cu impurity are expected to be occupied in the 



ground state, and the electrical neutrality of Cu impurity 
should be ensured by capturing two holes on Cu-related 
acceptor levels close to the top of the valence band, so 
that the resulting electron configuration can be denoted 
as Cu(d 10 p 2 ). Indeed such acceptor states were found in 
CaP:Cu samples, 12 although at that time the nature of 
these states remained unclear. 

Recently, ferromagnetism with a high Curie tempera- 
ture in p-type Cu-doped GaP was detected^ The EPR 
signal of the Cu 2+ state indicates that the 3d shells of Cu 
impurities are unfilled in this material in contradiction to 
the results of previous numerical calculations. This dis- 
crepancy gives us a motivation to revisit the problem. 

We present in this paper the results of numerical cal- 
culations of the electronic structure of Cu-doped GaP. 
Two different computation schemes are used, which give 
mutually complementary information about the behavior 
of weakly and strongly doped materials. The first one 
is the conventional local density approximation (LDA) 
scheme applied to the lattice of Cu^Gai-^P supercells. 
Similar methods were used for Mn^Gai-^Pn materials 
with Pn=As,N,Pi 14 i 15 The second method is based on 
the local Green function approach^ In this method the 
hybridization between the local impurity <i-orbitals and 
Bloch waves in the host semiconductor is calculated ex- 
actly, without any kind of artificial periodic boundary 
conditions, and approximations are made only when tak- 
ing into account the short-range part of substitution im- 
purity potential. 



II. GREEN FUNCTION APPROACH FOR 
ISOLATED IMPURITY 

A Green function calculation procedure based on 
the microscopic Anderson modelii was proposed three 
decades ag o 18 > 19 and later on summarized in Ref. 7. This 
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procedure deals with the local Green function 

G imp (r,r',z) = ^|A)(A|(z-H)- 1 |A)(A| . (1) 

A 

The set | A) includes both the electron states <^>f (r) of the 
electrons localized in the d-shell of impurity atom and the 
states ipb,-yfj,a (r) , which stand for "the Bloch tail" of the 
impurity wave function. These states describe the distor- 
tion inserted by a substitution impurity in the spectrum 
of a host crystal. They are superpositions of the Bloch 
waves, i/>b nt i<T = J2kn Ckn^kno-, where k an d n are the 
wave vector and the band index respectively, a is the 
spin quantum number. Here 7 is the index of the irre- 
ducible representation of the point group characterizing 
the symmetry of impurity and its surrounding, and \x de- 
notes its row. Therefore the function G; mp is diagonal 
in 7/i representation. The full Hamiltonian H includes 
the kinetic and potential energies of all electrons in the 
impurity atom and in the host crystal, as well as the 
Coulomb and exchange interactions between these elec- 
trons. The projection procedure (JTJ is exact in principle, 
and the poles of the Green function Gi mp describe both 
continuous and localized impurity related states in the 
doped crystal. In the practical realization of this method 
some approximations are unavoidable. The main simpli- 
fication, which we use here is the approximate treatment 
of the substitution potential 

AV(r - R ) = y cff (r - R ) - Vg(r - R ), 

where V^?(r — Ro) is the potential landscape for an 
electron in the host gallium atom in the site Ro and 
V e s(r — Ro) is the self-consistent potential for the elec- 
trons in the 3d shell of the Cu ion substituting for Ga in 
this site (see Section III for detailed definition of these 
potentials). We suppose that this potential is localized 
within the defect cell of the doped crystal. The "local 
substitution potential" approximation influences only the 
description of p-type acceptor states in the lower part of 
the forbidden energy band. It ignores possible contribu- 
tion of the Coulomb component of substitution impurity 
potential. This contribution is known to be small in the 
case of (Mn,Ga)As?i£ and one may hope for a similar sit- 
uation in (Cu,Ga)P. The principal advantage of the local 
substitution potential is that in this case the system of 
Dyson equations for the impurity-related components of 
the Green function (JTJ) defined as 

G w (z) = ( 7M |(z-iJ)- 1 | 7M ) (2) 

may be solved analytically 16 . It yields the equation 

G~l l {z) = z-E dl -M 1 {z)/Q{z), (3) 

for the d-electron Green function. The positions of elec- 
tron d-levels Sd-y are found self-consistently as a solution 
of the Schrodinger equation for Cu-related orbitals in the 
crystalline environment. The self energy in the right- 
hand side of Eq. contains two contributions. The 



term A4 7 (z) describes the hybridization between the d- 
orbitals and the band electrons 

a^HEt^ (4) 

nk 

where the hybridization integral is 

M 7 ,„ k = J ^ w (r)AV(r)^ nk (r)dr. (5) 

The energy bands e n k and Bloch functions 7/vik(r) of the 
host GaP crystal are calculated by means of the first 
principle full potential LAPW metho d 20 ' 21 (see Section 
III for details). 
The factor 

Q(z) = l-AV G° h (z). (6) 

in Eq. §5§ describes the short-range potential scattering, 
where 

AV ° = E / <kAV(r - Ro)^nk'dr (7) 

nknk' 

is the substitution impurity potential localized in the de- 
fect shell, 

G° h (z) = $>k|(* - ffo)- Vk) = £ (8) 

, , * c-nk 

nk nk 

is the single-site lattice Green function for the electrons in 
the non-doped host crystal described by the Hamiltonian 

As was shown in Ref. |Ty, the Green function §5§ de- 
scribes the hybridization between the impurity d-electron 
orbitals and the electrons in the imperfect host crystal, 
where the band electrons are influenced by the potential 
scattering AV. If this scattering is strong enough, it re- 
sults in splitting off of localized levels from the top of the 
valence band. This effect is also taken into account in ([3]): 
the positions of the corresponding levels before the hy- 
bridization are determined by zeros of the function Q(z) 
in the energy gap of the host crystal. 

One of the fundamental statements of the theory of 
transition metal impurities in semiconductors^ is the 
necessity to discriminate between the impurity levels in 
the gap obtained as solutions of a self-consistent mean- 
field Schrodinger equation for a doped crystal and the 
true addition/extraction energy of a d-electron to/from 
the valence/conduction band. The latter energies are de- 
termined by the energy balance of " Allen reactions" 6 i 7 i 22 

e n+1/n =e c - E{d n+1 ) + E(d n ) (9) 

Here E(d p ) is the total energy of doped crystal with the 
impurity having p electrons in 3d shell. Two Allen reac- 
tions describe the electron transition from the top of the 
valence band e v to the empty neutral (acceptor) level and 
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the electron transition from an occupied charged (donor) 
level to the bottom of the conduction band e c , so that the 
energies (O characterize the true positions of the impu- 
rity levels with respect to the band edges in the presence 
of strong Coulomb and exchange interactions. These en- 
ergies do not necessarily coincide with the mean-field so- 
lutions of the Schrodinger equation due to the violation 
of Koopmans' theorem for the impurity ions. 

To minimize the mismatch between the single-electron 
and many-electron states, Slater proposed a concept of 
"transition state". According to his arguments, the ion- 
ization energy for a state with n electrons in the 3d shell, 
E(d n ) — E{d n ~ 1 ) may be approximated by the energy 
e(n — 1/2) calculated within the LDA single-electron cal- 
culation scheme. More refined LDA+U approac h 23 ! 24 
takes non-Koopmans' corrections to the single-electron 
energies into account explicitly (although still semi- 
phenomenologically) . In terms of the Allen energies (O, 
the energy U is just the difference between and 
£«/«-!. We tegt below both thc lda+U method of 

Green function calculations and the standard LDA su- 
percell description of dilute (Gu,Ga)P semiconductor. 



III. IMPURITY GREEN FUNCTIONS IN 
LDA+U APPROXIMATION 

To realize a numerical version of the Green function 
method we use the local density approximation (LDA) 
and its modification LDA+U, which takes into account 
strong electron-electron correlations on the impurity site. 
This section outlines the application of the LDA+U 
method to systems with local defects with a particular 
emphasis on the transition metal impurities for which thc 
resonant scattering in the d (I = 2) channel plays a cru- 
cial role. Here we present only the principal features of 
the scheme leaving many more mathematical details and 
definitions in Appendix. In this section we retain the 
spin index, having in mind to use the spin-unrestricted 
LSDA+U version of this method for the calculation of 
spin properties of dilute magnetic semiconductors, al- 
though in the practical calculations below only the spin- 
independent LDA+U version is used. 



The LDA + U method incorporates a correction to the 
LDA energy functional which provides an improved de- 
scription of the electron correlations. The principal idea 
of the LDA + U method is to separate the electron sys- 
tem into two subsystems of the localized d-electrons for 
which the Coulomb interaction is accounted for by the 
Hubbard repulsion term hU^2 m > m , p m Pm' m the model 
Hamiltonian whereas the delocalized s- and p-electrons 
are described by an orbital independent one electron po- 
tential V lda {t). 

As a result the impurity Green function |T]) is defined 
by the Dyson equation 



where 



Gf{z) = G^'{z) l + Mf(z)Gf(*) 



MUz) = Ml(z)/Q(z). 



(10) 



(11) 



The Green function is a solution of this equation. We 
work in the spherically symmetric local basis (i = plm) 
instead of cubic harmonics expansion (i = pjfi) used in 
©. Here p is the index of repeating irreducible repre- 
sentations 7/i or Im, the analog of the principle quantum 
number n in a spherical atom. 

The bare d-elcctron Green function 



G f>{z) 



1 



z-ei- AU 



cr LDA+U 



(12) 



includes intraatomic correlations in the form of LDA+U 
potential consisting of three terms, 



Va, 



LDA+U 
ii.a 



rdc 



AV+^ + AV.v +AV a . 

ii,a 1 it, a ' it, a 



(13) 



Here the first term is the substitution LDA potential 



AU 



LDA 
m,plm-(7 



E G '»U"o / drr 2 AV^(r)R 2 plrT (r) , (14) 



The second term is the electron-electron interaction po- 
tential in the 3d shell, 



mmm m" 



u, 



)P~v 



mm"m"mj l J p l m " ,plm 



+ U mm m"m" Pplm" ,plm" 



(15) 



and the last term is the double counting compensation 
potential, parametrized as 



- ~UQ2 n plm,c ~\) + Aj2 n P^ ~\)- ( 16 ) 
mcr m 



Here we introduced the occupational matrix 

1 f eF 
P™t> = -~ lm / [ G W]w dz 

as a contour integral of the relevant matrix elements of 
the LDA+U Green function ([T2]). The Slater integrals 2 ^ 



4 




FIG. 1: (Color online) The embedded sphere and coordinate 
system used in our calculations. 



in the atomic limit read 

U mim2m3 m 4 = (mi,m 3 \V ee \m2,mi) 
21 

fe=0 

where the coefficients 



(17) 



a 



mim2m3m4 



4?T 



2k 



- ^2 (lmi\Y kn \lm2){lma\Y£ ri \lm4) 



— k 



can be expressed in terms of the Gaunt coefficients 

yl m+m 
T lrnl' m' 



The hybridization matrix elements (|5|) in the numera- 
tor of the mass operator now take the form 



M 



(see Appendix A.3). 



/ tf*(r)&V(r)^ PW (r-T s )e(r-T s )dr. (18) 



Here the Bloch wave functions PW are calculated by 
means of the linearized augmented plane wave (LAPW) 
method, t s is the vector connecting substitution impu- 
rity site taken as the point of origin with its nearest P 
neighbors in the zinc-blend lattice. 

In the impurity version of FLAPW method the defect 
site occupied by a Cu ion is surrounded by the "em- 
bedded sphere" with the radius r em b which includes the 
impurity sphere with the radius rjj (muffin-tin region, 
where the impurity potential is non-zero). Muffin-tin 
spheres r s = (r em b — r z?)/2 with a non-zero host lat- 
tice potential surround also the neighboring Ga sites (see 
Fig. [I]). The impurity centered basis set is chosen as 
a set of the linear augmented spherical wave (LASW) 
functions [see Eqs. (|A3|) . (|A4[) ]. In accordance with the 
LASW method, a set of Bessel functions is used in the 
remaining part of the sphere r om b. The wave functions in 
the two regions are matched by the standard boundary 
conditions imposed on the wave function and its deriva- 
tive. The Bloch functions f££ APW ( r ) of the host GaP 
crystal outside the embedded sphere are obtained by the 
self-consistent FLAPW method. Using the impurity cen- 
tered local LASW functions we calculate the matrix ele- 
ments of the host Green function projected onto the local 
spin polarized LAPW functions in the spherical intersti- 
tial site. After matching the boundary conditions (see 
Appendix), the matrix element (|18[) is transformed into 



J 



A 2 N 



MXi = ^^2v n (k e ) J2 T, ilY *m(K)G l Z-M",L"M" i drr 2 R; l (r)AV LflM 4r)^(k e ,r) LAPW (19) 



0=1 L" M" Ir. 



(g stands for the vectors of reciprocal lattice, see Ap- 
pendix). As was mentioned above, substituting Ga for a 
Cu impurity results also in an appearance of a potential 
component of the impurity potential, which is taken into 
account approximately by adopting the Koster-Slater- 
like single site scattering approximation^ Then in ac- 
cordance with Eq. (|10p . one may introduce the modified 
mass operator (fTTj) . where the zeros of the operator 
Q(z) © determine the impurity states, which arise in 



the doped crystal due to the potential scattering only. 

The scattering amplitude AVb is calculated by substi- 
tuting the LAPW wave functions ^ n 'k' LAPW (r) for the 
Bloch functions in Eq. Q . 

As a result the equation for the deep level energy de- 
termined as a pole of the impurity Green function ([3]) 
within the framework of the LDA+U technique reads 

z-e i -AV^ DA+u = MU^ (20) 
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It takes into account the resonance part of the scattering 
amplitude in the d (I = 2) channel and its mixing with 
the potential scattering states arising in the p (I = 1) 
channel^ 

The adspace augmentation^ is used to represent the 
Green function (or resolvent) G(z) for the GaP crystal 
with a Cu impurity in the matrix form, Eq. (lAljl . The 
impurity augmented Green function is subdivided into 
two blocks, of which the upper left corner block G° A (z) is 
constructed using the basis of i orbitals where i refer to 
the i-th state with the energy Si of the isolated adatom. 
The host is represented by the lower right corner block 
G° h (z). 

It is worth emphasizing that such a direct introduc- 
tion of the new adatom related states is very effective in 
the matrix formulation. Since the high energy part of 
the spectrum of the differential operator is well suited 
for the description of the strongly localized d-type im- 
purity states^ the issue of the necessary number of the 
host crystal bands becomes crucially important. The di- 
rect introduction of the d-states drastically simplifies the 
problem. The Dyson equation may be then split into two 
independently solvable equations (see Appendix) which 
finally allows one to carry out the calculations of the GaP 
host Green G°(z) using only 15 bands. 

The problem is treated self consistently, starting with 
the trial set of LAPW functions obtained with the help of 
the impurity potential, which in the zero's approximation 
is just a sum of the atomic potentials of the defect crys- 
tal. The self-consistency procedure for AT^(r) is carried 
out in a mixed fashion. The first two iterations use the 
arithmetic average scheme, which later on is effectively 
substituted by the Aitken schemed Just seven iterations 
produce the « 2 • 1CP 4 Ry self-consistency. 

The equations presented in this Section will be our 
working formulas for the LDA+U calculations of the 
Ga(Cu)P compound where the Cu atoms substituting 
Ga host atoms will be considered as isolated impurities. 
A possible exchange interaction between the Cu atoms 
and the resulting magnetic effects will be considered else- 
where. 



IV. DISCUSSION OF THE RESULTS 

This section presents the results of calculation of the 
electronic structure of Cu^Gai^^P obtained by means 
of the two methods, both using the LDA approxima- 
tion. The Green function approach is based on the band 
structure calculated by means of the FLAPW method 
discussed in the previous section. The supercell ap- 
proach uses the AS-LMTO method^ for the band calcu- 
lations. The Voskc 30 and Perdew-Wang 31 parametriza- 
tion scheme is used for the calculation of the exchange- 
correlation potential in the former and latter approaches, 
respectively. Brillouin zone (BZ) integration is performed 
using the improved tetrahedron method^ 

According to the present FLAPW and LMTO calcu- 
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FIG. 2: (Color online)The functions Re{G° h ) and lm{G° h ) 
GaP. 



fol- 



iations, the undoped GaP is a semiconductor with the 
1.83 eV FLAPW indirect gap and 1.61 eV ASA-LMTO 
gap between the top of the valence band (VB) at the T 
point and the bottom of the conduction band (CB) at 
the (0,0,0.875) point close to the X point of the fee BZ. 
A direct gap of 1.77 eV opens at the T point. The 6.8 eV 
width valence band is formed by the strongly hybridized 
P p and Ga s and p states while the states at the top of 
VB in the vicinity of T are formed by the P and Ga p 
states with a dominant contribution of the former. The 
band originating from the P related s states hybridized 
with the Ga related s states is found between —12.5 and 
—9.5 eV and separated by a gap of 2.8 eV from the bot- 
tom of the valence band. The density of states (DOS) of 
GaP is visualized in Fig. [5] as the imaginary part of the 
Green function G° ([8]) calculated by the LAPW method. 

A similar picture is obtained by direct band structure 
calculations within the ASA-LMTO method. The differ- 
ence in the widths of the energy gaps only weakly influ- 
ences the structure of the Cu-related states in the energy 
spectrum of doped samples. We start the discussion of 
these states with a discussion of the supercell calcula- 
tions. 



A. Supercell energy spectrum of Cui_ x Ga x P 

The electronic structure of Cu^Gai-^P with x varying 
from 0.125 down to ^0.016 was calculated using 2a x 
2a x 2a, 3a x 3a x 3a, and 4a x 4a x 4a supercells of the 
cubic zinc-blend lattice. Calculations for £=0.125 (1/8), 
0.063 (1/16), and 0.031 (1/32) were performed for Fi3m 
(216) fee, I43to bec (217), and P43m (215) simple cubic 
unit cells, respectively. The face-centered cubic cells with 
a = 3a and a = 4a allowed to simulate compositions 
with x w 0.037 (1/27) and x w 0.016 (1/64). In all 
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FIG. 3: (Color online) Bands calculated along some high 
symmetry directions and the total DOS for CuxGai-^P with 
a;=0.031. 

the calculations the Ga ion in the (0,0,0) position was 
substituted by the Cu ion with the same atomic sphere 
radius. This way the tetrahedral (T^) symmetry of the 
Cu impurity site was preserved. The positions of host 
atoms around the Cu impurity were not relaxed. 

Upon the Cu substitution Cu^Gai-^P becomes a 
metal with each Cu impurity creating 2 holes in the va- 
lence band. At all the compositions x studied in the 
present work the Fermi level (ef) crosses the three bands 
which are triply degenerate at the highest energy in the 
r point. At x=0.063 the top of the valence band lies 
0.42 eV above Ef and moves to 0.13 eV as the Cu 
concentration decreases to x=0.016. As an example, 
bands calculated along some high symmetry directions 
for Cu^Gai-ajP with x=0.031 are show in Fig. [3] At 
this Cu concentration the top of the valence band is sit- 
uated 0.22 eV above ef- 

Figure |H (lower panel) shows the density of Cu d states 
in CuzGai-zP with x = 0.031 projected onto the ir- 
reducible representations e and t2 of the T4 symmetry 
group. The densities of p states of the nearest (Pi and 
Gai) and next nearest (P2 and Ga2) neighbors of the Cu 
impurity are presented in the middle and upper panels 
of Fig. |1 

The calculations show that the Cu d shell is almost 
completely filled and the Cu valency is close to 1+. Cu d 
states of e symmetry (3z 2 — 1 and x 2 — y 2 ) form a DOS 
peak centered at —2.5 eV. They are completely occupied 
and do not contribute to the bands crossing the Fermi 
level. The main peak of the density of the £2 (xy, yz, 
and zx) states is located at —3 eV. However, another two 
peaks of t% DOS are clearly seen just at Ep and 0.5 eV 
below it. The origin of these peaks becomes more clear 
when the Cu ti DOS is compared to the density of p 
states of the Pi ion closest to Cu. The latter shows two 
prominent peaks exactly at the same energies. Similar 
peaks can also be observed in Gai DOS as well as in DOS 
of the more distant P and Ga ions not shown in Fig.[U An 
analysis of the partial occupations shows that of 2 holes 
(h) created by the Cu impurity only 0.18 h is provided 
by the Cu ti states. Another 0.48 h is distributed over 
the p states of 4 Pi ions whereas the remaining 1.36 h is 
spread over more distant neighbors. 
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FIG. 4: (Color online) A symmetry resolved density of Cu 
d states (lower panel) and the density of P p (middle panel) 
and Gap (upper panel) states calculated for Cu^Gai-^P with 
21=0.031. 



It is worth noting that in spite of the appearance of 
the narrow DOS peak exactly at ef the spin-polarized 
calculations failed to produce a ferromagnetic solution 
even for the highest Cu concentrations studied. Appar- 
ently, this can be explained by the delocalized character 
of the states responsible for the peak and an insufficient 
strength of the Hund's exchange coupling for P and Ga p 
states, which give the dominant contribution to the cor- 
responding bands. At the same time, the contribution of 
Cu d states, for which a strong on-site exchange interac- 
tion is expected, to the peak at Ef is relatively small. 

We also performed test calculations for a few values of 
x in E x G&i- x P, in which a Ga ion was substituted by 
a vacancy E. A vacancy creates one more hole in the 
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valence band as compared to Cu. Nevertheless, in the 
vicinity of the Fermi level the band structures calculated 
for E x Ga,i- x P are similar to those for Cu-doped GaP. 
In particular, the density of Pi states at and just below 
ef has the same two-peak shape. These peaks are also 
reflected in the density of E d states of the ti symmetry, 
however, they are much less pronounced than the corre- 
sponding peaks of Cu t 2 DOS. Significantly higher peaks 
can be observed in the density of E p states which also 
transform according to the t 2 representation. 

Thus, we may conclude that the bands crossing ef in 
Cu-rGax^P are mainly formed by the p states of the 
nearest to the Cu impurity Pi ions that split off from the 
top of the GaP valence band as a result of breaking of 
the covalent P p - Ga bonds at the impurity site. These 
states have t 2 symmetry and hybridize strongly with the 
corresponding Cu d states. These states are, however, 
rather extended, which leads to a relatively strong dis- 
persion of the split-off bands even for £=0.016. 
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FIG. 5: (Color online)Real and imaginary parts of mass op- 
erator Mi a (e) for (Ga,Cu)P. 



B. Cu-related energy states of isolated impurity 

Before turning to the calculation of the Cu-impurity 
related levels in the host GaP, let us look at the energy 
dependence of the self energy part (|4]) , which is respon- 
sible for the renormalization of 3d levels due to the hy- 
bridization with the host band states. The hybridization 
matrix elements M n k,i a are calculated by means of Eq. 
(fT8|) using the Cu 3c? impurity wave functions and LAPW 
functions of the GaP host. Since the LAPW wave func- 
tions are defined within the volume subdivided into two 
muffin-tin parts and the surrounding volume, the integra- 
tion in Eq. (fl8)) is carried out in all three parts separately 
accounting for all the hybridization contributions as well 
as for the covalency induced non-spherical components 
of the difference potential. 

Figure [5] represents the real and imaginary parts of 
Aii(s) obtained for the (Ga,Cu)P compound. Here in- 
dex i represents one of the components of the t 2 irre- 
ducible representation. Comparison of lmMt 2 (s) with 
the density of band states which is shown as Im(G^) 
in Fig. [2] demonstrates that the weighting of the den- 
sity of states with the squared hybridization matrix el- 
ement reproduces the general shape and van Hove sin- 
gularities of the partial p-component of DOS. The dif- 
ferences between Re(G°) and KeA4t 2 ( £ ) are more no- 
ticeable. Both these functions are sums of the Hilbert 
transforms of the DOS and weighted DOS for all the 
valence and conduction bands, respectively. Therefore 
these function not only map the singularities of DOS in 
the given band on the singularities of its Hilbert trans- 
form but also accumulate asymptotic contributions of 
higher and lower bands at the given e. This accumulation 
results in a noticeable smoothing of the A4t 2 ( e ) function 
in the —6 to eV range. Besides, weighting with M t 2 2 (e) 
strongly reduces the amplitude of ReA^t 2 (e) in compari- 
son with Re(G°). Such strong reduction means that the 



hybridization-induced renormalization of the atomic 3d 
levels of the isolated Cu impurity is small enough, and 
their positions are predetermined mainly by the impu- 
rity core potential and Coulomb interaction within the 
muffin-tin sphere ro- 

To compare the energy spectrum of the Cu impurity in 
GaP obtained by the Green function method with that 
given by LDA in the supercell calculation scheme, we 
first compute this spectrum by solving Eq. (|20[) within 
the LDA scheme without the second term AV,^ in the im- 
purity potential (|13p . Both the resonant and short range 
potential components of impurity scattering were taken 
into account. These calculations yield the value £ v — 0.66 
eV for the impurity dt 2 resonance in the valence band, 
which is higher than that in the supercell calculation, 
and the d e peak lies slightly above this level. Apparently, 
these peaks are related to the van Hove singularities in 
the heavy hole band. These resonances are shallower 
than those seen in the supercell DOS (Fig. [4]). As was 
mentioned above, the d e peak in the latter structure is 
located at e v — 2.5 eV. However, one should remember 
that the center of gravity of the valence band DOS is 
shifted downward with respect to its position in the pure 
GaP due to the transformation of d e and dt 2 levels into 
rf-bands (see Fig. [3]). Potential scattering built in the self 
energy part Mi(z) in Eq. (j2"0|) results in the appearance 
of an empty impurity level at e v + 0.168 eV. This accep- 
tor level may be identified with the x — ► limit for the P 
related p-structure at the top of the valence band in the 
supercell DOS (Fig. |H middle panel). The occupation of 
the impurity d-shell in this case is close to 10, like in the 
supercell calculations. 

The computation of the impurity spectrum within the 
LDA+U scheme yields a self-consistent solution for the 
electron spectrum only for the transition state 3d 8 5 of 
Cu impurity. This solution is described below. First, we 
determined the position of non-perturbed 3d-level of the 
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FIG. 6: Electronic structure of (Ga,Cu)P, calculated from Eq. (|20fl . The lowest of the five levels in the left panel correspond to 
the states (2, ±1), the next (l,m) levels are classified as (2,-2), (2,0) and (2,-1-2) (bottom-up). See the text for further discussion 



Cu atom and the correlation parameters U — J. The iso- 
lated impurity energy £j(+8.5) = —20.9 eV is calculated 
by means of the semi-relativistic RATOM program 33 
for the 3d 8 5 configuration, which corresponds to the 
concept of the transition state adopted in this paper. 
The intraatomic Coulomb repulsion of the d-electrons is 
treated in the LDA + U approximation and m-dependent 
Coulomb integrals (fT7|) are calculated. The choice of the 



Figure [6] depicts the electronic structure of (Ga,Cu)P 
calculated by the Green function method. We present 
here three versions of the calculations which account for: 
(a) resonant scattering, (b) short range potential scatter- 
ing, and (c) combined case. 

In the resonant scattering approximation (Fig. [5^) 
where the term Q~ 1 (z) is omitted in Eq. (|20|). there are 
four occupied levels in the valence band and one empty 
level in the energy gap. The occupied states correspond 
to the configuration Cu(d 8 ) of the impurity ion. These 
levels reflect the multiplet structure of this configuration. 
Although we used the orbital quantum numbers in our 
computation procedure, the calculated electron density 
distribution reveals the Td point symmetry of the impu- 



parameters U — 4.5 eV and J — 0.7 eV is based on 
the analysis of the occupation numbers in the transition 
state approach^ The self-consistent single electron 3d- 
level for the embedded Cu impurity in the 3d 8 5 configu- 
ration is in resonance with the valence band of GaP host 
crystal, and the impurity-related resonance and discrete 
states are found as solutions of Eq. (12"0"1) . 



I 

rity surrounding. In terms of the corresponding cubic 
harmonics the lowest state has the ti symmetry, the two 
next levels belong to the e-representation, and the empty 
state in the energy gap is the t 2 state of the configuration 
d 9 . In terms of the Allen diagrams ((9]) these levels cor- 
respond to the addition energy e\ 8 = -E(e 4 i|) — -E(e 4 if ) 

and £e^ 8 = E(e 4 t^) — E{e 3 t\) for d t2 and d e quan- 
tum numbers, respectively (see similar classification for 
(Ga,Ni)P in Ref. [ill). The final 3d 8 states belong to the 
3 T 2 (F) and 3 Ti(F) representations in the Tanabe-Sugano 
classification£L±L±£ 

The energy interval between the multiplet of occupied 
levels in the valence band and the empty level in the en- 
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ergy gap is < 4eV, which is comparable with the value 
of the input parameter U — J = 3.8 eV. The hybridiza- 
tion renormalization due to the self energy M. (z) in Eq. 
plj) is 0.115 eV for the occupied levels and 0.182 eV for 
the empty level. In the calculation procedure described 
above, the difference in hybridization shifts for t% and 
e-levels was neglected, because the hybridization (ligand 
field) contribution is small enough for the Cu impurity 
ion. 

Figure [6Jd exhibits the net contribution of poten- 
tial scattering ([TJJ to the formation of impurity-related 
states. The levels shown in this figure are obtained from 
PPJl with M substituted for Q' 1 [see Eq. (TJT])]. The 
states in the occupied part of the spectrum are the im- 
purity resonances in the valence bands around the max- 
ima of the partial p-wave contributions at the energies 

6eV and 2eV (cf. Figs. Hand E]). The p-level 

arises at the energy +0.68eV above the top of the valence 
band. 

Both the d- and p-like states are found in the solution 
of Eq. (HOD with the full self energy M (Fig. [3b). The 
most significant difference between the combined spec- 
trum of Fig. [5b and those of Fig. [5ji,b is the noticeable 
hybridization between p and d t2 resonances in the valence 
band, whereas thee d e levels are only slightly shifted. The 
shallow p-level in the energy gap is pinned to its original 
position shown in Fig. [5b>, in spite of dp hybridization. 
All these results agree with qualitative predictions of the 
analytical model taking into account both resonant and 
short-range potential stattering^ 

There is no straightforward way to compare the results 
of LDA+U calculations with those obtained within the 
LDA scheme, since the former method uses the fitting 
parameters U, J, whereas the latter one is based on the 
variational approach, which formally gives the solution 
corresponding to the minimal total energy. We only may 
estimate the total energies of the two solutions by com- 
paring the positions of the impurity levels obtained by 
both methods within the same Green function approach. 
LDA procedure gives the occupied e and ti levels at the 
energies ~ e„ — 0.64 to 0.66 eV below the top of the va- 
lence band and the shallow p-level at the energy e^+0.168 
eV, which corresponds to the configuration d 1Q p 2 : two 
holes neutralize the excess charge in the d shell, which 
means that the triply degenerate p-level is occupied by 
one electron. In the LDA+U solution the occupied t-i and 
e-levels lie essentially deeper in the valence band at the 
energies ~ e v — 2.3 to 1.8 eV, Cu ion behaves as the iso- 
electronic impurity Cu 3+ (d 8 ), and the acceptor p-levels 
are triply occupied in the neutral impurity state. The 
comparison of single-electron energies for the two solu- 
tions gives the energy gain ~ 10.7 eV for the latter state. 
It is hardly probable that the exchange-correlation con- 
tribution may change the energy balance in favor of a 
state with the fully occupied 3d shell of the Cu impurity. 

Comparing the electronic structures of (Ga,Cu)P ob- 
tained by the supercell and Green function methods, one 
may indicate both similarities and dissimilarities in the 



description of impurity-related states. 

First, both methods provide the same mechanism for 
formation of the shallow p-levels in the energy gap of 
the host material, which merge into the impurity band 
at a high enough dopant concentration. These levels are 
split off from the top of the valence band and partially 
hybridized with the ti levels in the valence band. 

Second, the spectral density of the impurity related 
<i 7 -states is concentrated mainly in the valence band 
with the e?t 2 component lying below the d e component. 
Here, however, the important difference between the two 
approaches should be emphasized. As was mentioned 
above, the ei 7 resonances calculated within the Green 
function LDA approximation are shallower than those 
found by means of the supercell approach. One may see 
also a difference in the d e — d t2 splitting: it can be es- 
timated as ~ 0.5eV in the supercell calculations and as 
< 0.3eV in the Green function calculations. The main 
reason of this difference is the fact that the impurity 3d- 
levels are transformed in an effective rf-bands in the peri- 
odic supercell structure, and the hybridization repulsion 
between the two Bloch waves is stronger than that be- 
tween the localized d-levels and periodic partial p-waves 
in the Green function approach. The same argument is 
valid for the quasiband method used in the calculations 
of Ref. [n]; where the Cu-related d - levels are located 
even deeper than in our supercell calculations at the en- 
ergies ~ 3 — 4.5 eV below the top of the valence band of 
GaP. 

The most important difference between the results de- 
scribed in Subsections IVA and IVB is of course the differ- 
ence in the electron configuration of Cu impurity, which 
is d w p 2 in the supercell calculations and d 8 in the Green 
function calculations. Available experimental data^ are 
in favor of the configuration d 9 p. At this stage we have 
no exhaustive explanation of these discrepancies. First, 
our scheme should be extended to the spin-unrestricted 
LSDA solution and to the multiimpurity case. We ex- 
pect that the charge configuration of Cu ions is highly 
sensitive both to the spin state and to the interimpu- 
rity coupling. Second, more experimental investigations 
are necessary, which would reveal the role of concomi- 
tant defects, the annealing conditions, the thickness of 
the film and other technological factors. It is also worth- 
while checking whether the use of LDA+U method in the 
supercell approach may result in the configuration with 
an incomplete 3d shell of the Cu impurity. We leave all 
these questions for further investigations. 



V. CONCLUDING REMARKS 

The numerical solution of the Dyson equation (j2"0|) de- 
rived by means of the Green function method reveals sim- 
ilarities and dissimilarities between the electronic struc- 
tures of the Mn impurity (half-filled 3d shell in atomic 
state) and Cu impurity (completely filled 3d shell in 
atomic state) substituting for Ga in zinc blende semi- 
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conductor. Our calculations show that unlike Mn, which 
retains its stable half-filled 3d 5 shell in the host GaAs and 
GaP crystals,— the Cu impurity may release some of its 
d-electrons from the stable filled shell 3d 10 to minimize 
the total energy of doped crystal, at least in the wide-gap 
GaP. Our theoretical result partially agrees with the ex- 
perimental observation of Cu ions with unfilled 3d shell 
in GaPfi£ It paves a way to theoretical explanation of the 
ferromagnetic ordering in Gai-aCu^P crystals, although 
for this purpose further development of the Green func- 
tion method is necessary. The results of the numerical 
study of magnetic ordering by means of the Green func- 
tion method will be published elsewhere. 

This work is partially supported by the Max-Planck 
Gesellschaft during the stay of O.F., K.K. and V.F. in 
MPIPKS (Dresden) , where this work was completed. 



APPENDIX A: DETAILS OF COMPUTATIONAL 
SCHEME 

In order to realize the GF approach in a compu- 
tational scheme we make use of the local density ap- 
proach (LDA)^ and its LDA+U modification^ which 
accounts for a strong electron-electron interaction. The 
approximation 3 ^ is used for the exchange-correlation po- 
tential. The band structure of the GaP semiconductor is 
calculated by means of the ab-initio full potential all elec- 
trons LAPW method. 21 This method presents the charge 
density and the crystal potential as a series of the spher- 
ical harmonics inside the muffin-tin spheres and of the 
plane waves outside the spheres. The self-consistent elec- 
tronic band structure is determined by solving a single 
particle Dirac equation by using the variational method 
in LAPW - function basis PW ( r )}- In order to eval- 
uate the Coulomb part of the crystal potential we use the 
concept of multipole potentials and solve the Dirichlet 
problem for the sphere with all the contributions being 
treated on equal footing^ The exchange-correlation po- 
tential is approximated by the Pade approximant tech- 
nique in order to interpolate accurately the recent Monte 
Carlo results with the RPA spin-dependent data^ The 
Fourier components of the exchange - correlation poten- 



tial in the interstitial region are fitted in the least square 
method by applying the singular value decomposition 
procedure. The charge density in the interstitial region 
is calculated in ca. 2000 to 3000 random points in the 
irreducible wedge of the Wigner-Seitz cell. 

In order to find the self energy Aii(z), one has to calcu- 
late the matrix elements M„k,i between the band states 
|nk) and the states \i) = \plm) of the impurity atom. 
A computational scheme based on the augmented Green 
functions 38 



G°(*) 



G° A (z) 



Gl(z) 
is developed for this sake. Here 



(Al) 



G>,r';z)= 



&(r)#(r) 



i— p,l,m 

is the impurity Green function, whereas the host crystal 
is represented by 



G° h (r,r';z) = Y: £ " k ( 1 " k ^ 



n=l WelBZ 



The wave functions of electrons localized in the im- 
purity 3d shell are defined within the impurity sphere 
r <r D (see Fig. <^ m (r) = R pl (r)Yi m (r). The radial 
parts of these functions are defined as solutions of the 
equation 

[-V 2 + V?(r) + AV(r)] R pl (r) = e pl R pl (r), (A2) 

and the angular parts are represented by the spherical 
harmonics. The Bloch wave functions are expanded in 
the reciprocal wave vectors k e = k + TC e 

N 

0=1 

where 



Vk e (r) 



'Ho 



with 



The following notations has been used above: 



a\ s) (k e )n\ s) (s h p) + b\ s \k e )n\ s) (s h p) Y lm (p) 



e int (r) 



1, re flint — volume of the 
interstitial region 

0, otherwise, 
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e s (p) = e s (r-r s ; 



1, p £ fl s — volume of the 
s sphere region 

0, otherwise, 



v n (k e ) are eigenvectors of LAP W variation procedure; n 
is the number of the accounted energy bands, flo is the 
volume of the Wiener- Seitz cell, a, (s) (k e ) and &| s) (k e ) 
are the muffin-tin coefficients in the LAPW method, 
and TZ^ (si , p) = ^Tl\ s \e, p)\ £l for the fixed energy £/. 
iz\ s ^ (ei , p) is the radial part of the LAPW function. 



1. Choice of the localized basis 

Calculations of the electronic structure of defects in 
crystals are usually based on the pseudopotential or 
LCAO + pseudopotential approac h 39 i 40 . This method 
requires a large number of the Gaussian orbitals and cal- 
culation of their overlap integrals. Instead we perform 
here an all-electron calculation which allows one to realize 
the spin-polarization LDA + U scheme. This approach 
uses the basic set of No functions 



Xm( f ) = XpLM{r,tf,<p) 



in order to obtain the orthonormal basis 
X,W=E( L_1 )L'^(r). 

Then the Green function of the host crystal is projected 
onto the localized basis 



A I 



n=l kelBZ 



(x»\e\nr w )(n£ PW \®\x 

z - £„k 



and calculated by means of the analytical tetrahedron 
method 32 within the irreducible part of the Brillouin zone 
(IBZ) 



n=l kElBZ 



Here 



Tn^v = (X»\®\<£ PW ){<£ PW \Q\X»)- 



The coefficients ui„k(-z) depend only on the dispersion 
relation e n k and can be computed only once. 



2. Self energies for impurity Green function 



F pL (r)Y LM (?) , for r <r D , 
jL(K p Lr)Y LM (r) for r > r D 



(A3) 



Here L is a non- negative quantum number and — L < 
M < L, the inverse length k p l is defined by zeros of the 
Bessel function jjj(KL p r em b) = for the radius r em b of the 
embedded sphere; p is the integer number enumerating 
these zeros. The radial part of the wave function (|A3[) is 



F pL (r) = a L (n pL )R L {e L ,r) + b L (n pL )R L (s L ,r). (A4) 
Here the parameters 



cll(k p l) 



are used to match the function (|A4[) to the Bessel func- 
tions outside the muffin-tin region, = R^{eL,r), 



RlR'l — RlR'l 



R' L jL{KpLrp) - R' L i' L {K p Lrp) 
RlR'l — RlR'l 



R'r 



dR L {e L ,r) 



dr 



Rl 



dR L (e L ,r) 



dej 



The above basis x^( r ) was used in the Cholesky de- 
composition S = L • for the overlap matrix 



xl{?)Xv{*)dv 



The impurity Green function (|10[) contains several self 
energy corrections to the atomic levels Si . Two of them 
given by Eqs. (p~5|) ■ (fl6l) arising from the Coulomb inter- 
action are responsible for the multiplet structure of the 
energy levels, potential contribution (TT3]) results in the 
crystal field splitting of these levels, and the resonance 
self energy (fTTj) is the analog of ligand field correction in 
conventional theory of transition metal impurities^ This 
section discusses the calculation of the two last terms 
within the Green function formalism. 

The resolvent operator AG(z) and the corresponding 
density variation Ap(r) is calculated both for the host 
block [AG(z), Ap(r)] and for impurity block [AGu(z), 
Api (r)] of the secular matrix (|A1[) . When calculating the 
contour integrals resulting in (|A7[) we use semi-circular 
contour from the bottom of the valence band Sb to the 
Fermi energy e^. The charge dependent difference po- 
tential AV(r) is not necessarily spherically symmetric. 
We define the substitution impurity potential as the dif- 
ference 



AV[p(r)] = V eff [p(r)]-V£[pUr)] 



(A5) 



between the true self-consistent effective potential 
V e ff[p(r)] and the effective self-consistent potential 
V^[p°(r)] of the host crystal, both taken in the LDA 
approximation. Here p(r) and p°(r) are the respective 
electron densities. 

The impurity correction to the host Green function of 
the crystal induced by the potential (Q2 



AG(z) = G(z) 
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is found from the corresponding Dyson equation 41 
AG(z) = 



l-G^-AV-O-- 1 )^ '-I 



Here 



X*(r)AV\p(r)] x Ar)dT 



and I is a unit matrix. 

The density variation is calculated using the equation 

N D N D 

Ap(r)=Im^^Ap^ X/i (r) X :(r) (A6) 

fl — l u — l 



where 



and 



A^=((L- 1 )tA /3 L- 1 )^ 



Ap = -- [ F AG(z)dz. (A7) 

The lower integration limit Sb is chosen to include all 
the relevant band and impurity states, ep is the Fermi 
energy. To compute the integral (|A7[) . we introduce the 



contour C in the complex plane z enclosing all the poles 
up of the Green function up to the Fermi energy in the 
charge density integration. 

With Ap(r) calculated by means of Eqs. (|A6[) and 
(|A7[) we calculate anew the charge dependent impurity 
potential 



AV(r) = kV LM {r)Y LM {v) 



(A8) 



LM 



in the " embedded cavity" , which is not spherically sym- 
metric. The density variation can be similarly repre- 
sented as 



Ap(v) = Y J ^PLM(r)Y LM (v) 



(A9) 



LM 



where one readily obtains 

Ap L "M"{r) 



J2Y. A PpL,P'L'r PP >(L,L';r) £ G%%+tf M „ (A10) 

pp> LL' M=-L 

with 



T ppl (L,L';r) = 

aL(K P L)aL'(Kp'L)RL(£L,r)R L ,(e L ,r) + a L {K pL )b L >{K pL )RL(e L ,r)Ry{e Ll r) + 
aL'(KpL)bL(n P L)RL(£L,r)R L >(e L >,r) + b L (K pL )b L >(K p > L >)R L (s L ,r)R L ,(s L ,r), for < r < r D 



, jh{npLr)jL'{np'L'r) 



for td < r < r, 



ernb 



(All) 



and 



G 



being the Gaunt coefficients. 

Next we separate the impurity and host parts in the 
density correction 



Ap L „ M „(r) = Ap L „ M „{r) + Ap { °), u „(r) 



where 



N D 



Ap L „ M „(r) = ^2 Ap^r pp (L, L;r)G 



LM 

LML"M" 



li=l 



(L, L'\ t)G\ ml „ M ii 

li=2 fj,' = l 



is the host contribution, and 



Ap^, M „(r)=R 2 pl (r) £ Ap lm , lm G 



lm+M" 
ImL" M" 



-—I 



is the substitution impurity contribution. The functions 
Rpi (t) are the radial parts of the impurity centered local 
orbitals (|Al|) 

Using this density correction, we calculate the impurity 
related self energy AV i f DA 



AV; 



LDA 



lm-M"L"M" 



L" M" 



drr 2 AV L „ M »(r)R 2 pl (r) 



(A12) 
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and substitute it into the Green function (jTSJ) . the integral (fl8|) and matching the boundary conditions 

Self energy correction At; in pip contains the matrix in accordance with the procedure described above, the 
elements (fT5| . After substituting the potential (|A5|) in hybridization matrix element acquires the form (fl9|) . 



1 T. Jungwirth, J. Sinova, J. Masek, J. Kucera, and A.H. 
MacDonald, Rev. Mod. Phys. 78, 809 (2006). 

2 R. Bouzerar, G. Bouzerar, and T. Ziman, Phys. Rev. B 
73, 024411 (2006). 

3 M. Strassburg, M.H. Kane, A. Asghar, Q. Song, Z.J. 
Zhang, J. Senawiratne, M. Alevi, N. Dietz, C.J. Summers, 
and I.T. Ferguson, J. Phys.: Condens. Matter 18, 2615 
(2006). 

4 J.M.D. Coey, M. Venkatesan, and C.B. Fitzgerald, Nature 
Materials 4, 173 (2005). 

5 K.A. Griffin, A.B. Pakhomov, CM. Wang, S.M. Heald, 
and K.M. Krishnan, J. Appl. Phys. 97, 10D320 (2005); 
Phys. Rev. Lett. 94, 157204 (2005). 

6 A. Zunger in Solid State Physics, edited by F. Seitz and 

D. Turnbull (Academic Press, New York, 1986), Vol. 39, 
p. 275. 

7 K.A. Kikoin and V.N. Fleurov, Transition Metal Impuri- 
ties in Semiconductors (World Scientific, Singapore, 1994). 

8 A.K. Bhattacharjee, C. Benoit a la Guillaume, Solid St. 
Commun. 113, 17 (2000) 

9 P.M. Krstajic, F.M. Peeters, V.A. Ivanov, V. Fleurov and 
K. Kikoin, Phys. Rev. B 70, 195215 (2004) 

10 T. Jungwirth, J. Sinova, A.H. MacDonald, B.L. Gallagher, 
V. Novak, K.W. Edmonds, A.W. Rushforth, R.P. Cam- 
pion, C.T. Foxon, L. Eaves, K. Olejmk, J. Masek, S.-R. 

E. Yang, J. Wunderlich, C. Gould, L.W. Molenkamp, T. 
Dietl, and H. Ohno, Phys. Rev. B 76, 125206 (2007) 

11 V.A. Singh and A. Zunger, Phys. Rev. B 31, 3729 (1985) 

12 L.-A. Ledebo and B.K. Ridley, J. Phys. C: Solid State 
Phys. 15, L961 (1982) 

13 A. Gupta, F.J. Owens, K.V. Rao, Z. Iqbal, J.M. Osorio 
Guille, and A. Ahuja, Phys. Rev. B 74, 224449 (2006); 

F. J. Owens, A. Gupta, K.V. Rao, Z. Iqbal, J.M. Osorio 
Guille, A. Ahuja, and J.-H. Guo, IEEE Trans. Magn. 43, 
3043 (2007). 

14 A.B. Shick, J. Kudrnovsky, V. Drchal, Phys. Rev. B 69, 
125207 (2004). 

15 L. Kronik, M. Jain, and J.R. Chelikowsky, Appl. Phys. 
Lett. 85, 2014 (2004). 

16 V.N. Fleurov and K.A. Kikoin, J. Phys. C 19, 887 (1986). 

17 P.W. Anderson, Phys. Rev. 124, 41 (1961). 

18 V.N. Fleurov and K.A. Kikoin, J. Phys. C: Solid State 
Phys. 9, 1673 (1976). 

19 F.D.M. Haldane and P.W. Anderson, Phys. Rev. B 13, 
2553 (1976). 

20 L.F. Mattheiss and D.R. Hamann, Phys. Rev. B 33, 823 



(1986). 

21 O.V. Farberovich, S.V. Vlasov, K.I. Portnoi, and A.Yu. 
Lozovoi, Physica B 182, 267 (1992). 

22 J.W. Allen in Proc. 7-th Int. Conf. Physics of Semiconduc- 
tors (Paris, Dunod, 1964), p. 781. 

23 V.I. Anisimov and J. Zaanen, O.K. Andersen, Phys. Rev. 
B 44, 943 (1991). 

24 I.V. Solovyev, P.H. Dederichs, V.I. Anisimov, Phys. Rev. 
B 50, 16861 (1994). 

25 M.T. Czyzyk and G.A. Sawatzky, Phys. Rev. B49, 14211 
(1994). 

26 A.R. Williams, P.J. Feibelman and N.D. Lang, Phys. Rev., 
B26, 5433 (1982). 

27 U. Lindefelt and A. Zunger.Phys. Rev., B26, 846 (1982). 

28 A.C. Aitken, Proc. Roy. Soc. Edinburgh, 46, 289 (1926). 

29 A.R. Williams, J.Kubler, CD. Gelatt, Phys. Rev., B19, 
6094 (1979). 

30 S.H. Vosko, L. Wilk, and M. Nussiar, Can. J. Phys., 58, 
1200 (1980). 

31 J.R Perdew and Y. Wang, Phys. Rev., B45, 13244 (1992). 

32 P. Lambin and J.R Vigneron, Phys. Rev., B29, 3430 
(1992); RE. Blochl, O. Jepsen, and O. K. Andersen, Phys. 
Rev., B49, 16223 (1995). 

33 O.V. Farberovich, S.V. Vlasov, and CP. Nizhnikova, 
Program of the self-consistent relativistic calculation of 
the atomic and ionic structures in LSDA approximation, 
VINITI No. 2953-83 (1983) (Russia). 

34 K. Sato, P.H. Dederichs, H. Katayama-Yoshida and J. Ku- 
drnovsky, J. Phys.: Condens. Matter., 16, S5491 (2004). 

35 C. Timm and A.H. MacDonald, Phys. Rev. B 71, 155206 
(2005) 

36 O.V. Farberovich, Electronic structure and physical prop- 
erties of compounds with d- and f-metals, Ph.D. thesis, 
Voronezh University (Russia), 1984. 

37 V.I. Anisimov, I.V. Solovyev, M.A. Korotin, M.T. Czyzyk, 
and G.A. Sawatzky, Phys. Rev. B48, 16929 (1993). 

38 H. Katayama-Yoshida and A. Zunger, Phys. Rev., B31, 
7877 (1985). 

39 J. Bernholc, N.O.Lipari and S.T. Pantelides, Phys. Rev. 
Lett., 41, 895 (1978). 

40 G.A. Baraff and M. Schluter, Phys. Rev. Lett., 41, 892 
(1978). 

41 G. Wachutka, A. Fleszar, F. Maca, and M. Scheffler, J. 
Phys.: Condens. Matter., 4, 2831 (1992). 



